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1. INTRODUCTION 

As part of the ongoing effort at NASA to improve both the durability and reliability of 
hot section Earth-to-Orbit engine components, continuing improvements must be made in 
existing finite element and finite difference methods, and alternate techniques, such as the 
boundary element method, must be explored. Despite this considerable effort, the accurate 
determination of transient thermal stresses in these hot section components remains one 
of the most difficult problems facing engine design/analysts. For these problems, the 
temperature distribution is strongly influenced by the external hot gas flow, the internal 
cooling system, and the structural deformation. Currently, experimentally-determined film 
coefficients and ambient temperatures are required for use as boundary conditions for the 
thermal stress analysis of the structural component. The determination of these coefficients 
is obviously an expensive and time-consuming task. Very recently an attempt was made 
by Gladden (1989) to use a finite difference-based Navier-Stokes code to approximate the 
thermal boundary conditions, and to then input these into a finite element structural 
analysis package. However, the most effective way to deal with this problem is to develop 
a completely integrated solid mechanics, fluid mechanics, and heat transfer approach. 

In the present work, the boundary element method (BEM) is chosen as the basic anal- 
ysis tool principally because the critical surface variables (i.e., temperature, flux, displace- 
ment, traction) can be very precisely determined with a boundary-based discretization 
scheme. The price that must be paid for this precision is that any BEM formulation re- 
quires a considerable amount of analytical work, which is typically absent in the other 
numerical methods. 

This report details the progress made, during the period November 1988 - November 
1989 in a multi-year program commencing in March 1986, toward the development of a 
boundary element formulation for the study of hot fluid-structure interaction in Earth-to- 
Orbit engine hot section component?. Most of the work reported in previous years under 
this program was directed toward the examination of fluid flow, since boundary element 
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methods for fluids are at a much less developed state. This year significant strides have 
been made, not only in the analysis of thermoviscous fluids, but also in the solution of the 
fluid-structure interaction problem. 

During the first half of this past year, the convective viscous integral formulation was 
derived and implemented in the general purpose computer program GP-BEST. The new 
convective kernel functions that were developed as a part of this effort, m turn, necessi- 
tated the development of refined integration techniques. As a result, however, since the 
physics of the problem is embedded in these kernels, boundary element solutions can now 
be obtained at very high Reynolds number. Flow around obstacles can be solved approxi- 
mately with an efficient linearized boundary-only analysis or more exactly by including all 
of the nonlinearities present in the neighborhood of the obstacle. 

The other major accomplishment of 1989 has been the development of a comprehensive 
fluid-structure interaction capability within GP-BEST. This new facility is implemented 
in a completely general manner, so that quite arbitrary geometry, material properties and 
boundary conditions may be specified. Thus, a single analysis code (GP-BEST) can be 
used to run structures-only problems, fluids-only problems, or the combined fluid-structure 
problem. In all three cases, steady or transient conditions can be selected, with or without 
thermal effects. Nonlinear analyses can be solved via direct iteration or by employing a 
modified Newton-Raphson approach. 

In the next section, a brief review of the recent applicable boundary element literature 
is presented. This is followed by the development of integral formulations for the ther- 
moelastic solid in Section 3 and for the thermoviscous fluid in Section 4. A number of 
detailed numerical examples are included at the end of these two sections to validate the 
formulations and to emphasize both the accuracy and generality of the implementation. 
Then, in Section 5, the fluid-structure interaction facility is discussed. Once again, several 
examples are provided to highlight this unique capability. Section 6 contains a collection 
of potential boundary element applications that have been uncovered as a result of work 


2 



related to the present grant. For most of those problems, satisfactory analysis techniques 
do not currently exist. The remaining sections summarize the progress achieved to date, 
and outline the work plan for the next year. Tables and figures appear at the end of each 
section, while references are provided in Appendix A. 
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2. LITERATURE REVIEW 

Virtually nothing has appeared in the literature on the analysis of coupled thermo- 
viscous fluid-structure problems via the boundary element method, except for Dargush 
and Banerjee (1988,1989a) which is a summary of early work performed under this grant. 
However, a number of publications have addressed the fluid and structure separately. 

In general, the solid portion of the problem has been addressed to a much greater 
degree. For example, a boundary-only steady-state thermoelastic formulation was initially 
presented by Cruse et al (1977) and Rizzo and Shippy (1977). Recently, the present 
authors developed and implemented the quasistatic counterpart (Dargush, 1987; Dargush 
and Banerjee, 1989b), which is presented in detail in Section 3. Others, notably Sharp 
and Crouch (1986) and Chaudouet (1987), introduce volume integrals, to represent the 
equivalent thermal body forces. A similar domain based approach was taken earlier by 
Banerjee and Butterfield (1981) in the context of the analogous geomechanical problem. 

An extensive review of the applications of integral formulations to viscous flow prob- 
lems was included in a previous annual report (Dargush et al, 1987), and will not be 
repeated here. Interestingly, only a few groups of researchers are actively pursuing the 
further development of boundary elements for the analysis of viscous fluids. The work re- 
ported in Piva and Morino (1987) and Piva et al (1987) focuses heavily on the development 
of fundamental solutions and integral formulations with little emphasis on implementation. 
On the other hand, Tosaka and Kakuda (1986, 1987), Tosaka and Onishi (1986) have im- 
plemented single region boundary element formulations using approximate incompressible 
fundamental solutions. This latter group has developed sophisticated non-linear solution 
algorithms, and consequently, are able to demonstrate moderately high Reynolds number 
solutions. 

The most recent work from the above researchers has been collected into a volume en- 
titled Developments in BEM - Volume 6: Nonlinear Pro blems of Fluid Dynamics, edited 
by Banerjee and Morino. Contributions from Wu and Wang, and Bush and Tanner are 
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also included, along with two chapters from the present co-authors. The volume, pub- 
lished by Elsevier Applied Science Publishers with availability in mid- 1990, will provide a 
state-of-the-art review of boundary element fluid dynamics. However, the convective ther- 
moviscous formulations of Section 4 are a significant further advancement which permit 
solutions for high Reynolds number flows. 


5 



3. INTEGRAL FORMULATION FOR SOLIDS 

3.1 Introduction 

In the current section, a surface only time domain boundary element method (BEM) 
will be described for a thermoelastic body under quasistatic loading. Thus, transient heat 
conduction is included, but inertial effects are ignored. This BEM was first developed 
as part of the work performed during the second year (1987) of this grant. Since that 
time a number of improvements and extensions have been incorporated. During 1989, 
the algorithms for numerical integration have been made more efficient as well as more 
accurate, and a comprehensive PATRAN interface has been added to aid in the post- 
processing of the boundary element results. Additionally, a streamlined approach for 
uncoupled thermoelasticity was introduced (Dargush and Banerjee, 1989b). 

Details of the integral formulation for 2D plane strain is presented below. Separate 
subsections present the governing differential equations, the integral equations, an overview 
of the numerical implementation, and a couple of simple examples. Similar formulations 
have also been developed for three-dimensional (Dargush and Banerjee, 1990a) and ax- 
isymmetric problems (Dargush, 1987). 


3.2 Governing Equations 

With the solid assumed to be a linear thermoelastic medium, the governing differential 
equations for transient thermoelasticity can be written 


(A + n) 


d 2 Uj 

dx^ d xy 


d 2 Ui 


dO 




d$ d 2 6 
P Ct dt dijdxj 


(3.1a) 


(3-16) 


where 

ui displacement vector 
6 temperature 
t time 

Xi Lagrangian coordinate 
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k thermal conductivity 


p mass density 

c« specific heat at constant deformation 
A, p Lame constants 

a coefficient of thermal expansion 

Standard indicial notation has been employed with summations indicated by repeated 
indices. For two-dimensional problems considered herein, the Latin indices t and j vary 

from one to two. 

Note that (3.1b) is the energy equation and that (3.1a) represents the momentum 
balance in terms of displacements and temperature. The theory portrayed by the above 
set of equations, formally labeled uncoupled quasistatic thermoelasticity, can be derived 
from thermodynamic principles. (See Boley and Weiner (1960) for details.) 


3.3 Integral Representations 

Utilizing equation (3.1) for the solid along with a generalized form of the reciprocal 
theorem, permits one to develop the following boundary integral equation: 


«*»(£)«*(£.*)■ f 90a * tp(X,t) - Ua*^0{X, t)^dS(X). 
J 8 


( 3 . 2 ) 


where 

a, p indices varying from 1 to 3 
s surface of solid 

u a ,t a generalized displacement and traction 

tia = [«1 “2 0] T 

t a = [ ( 1 *2 SF 

0 ,q temperature, heat flux 

g a 0,f a p generalized displacement and traction kernels (Dargush, 1987,1989b) 
c a 0 constants determined by the relative smoothness of s at £ 


7 



and, for example 


9a * ta — [ 9afi( x > £> r )^a( x i r )dT 

Jo 


denotes a Riemann convolution integral. 

In principle, at each instant of time progressing from time zero, this equation can be 
written at every point on the boundary. The collection of the resulting equations could then 
be solved simultaneously, producing exact values for all the unknown boundary quantities. 
In reality, of course, discretization is needed to limit this process to a finite number of 
equations and unknowns. Techniques useful for the discretization of (3.2) are the subject 
of the following section. 

3.4 Numerical Implementation 

3.4.1 Introduction 

The boundary integral equation (3.2), developed in the last section, is an exact state- 
ment. No approximations have been introduced other than those used to formulate the 
boundary value problem. However, in order to apply (3.2) for the solution of practical en- 
gineering problems, approximations are required in both time and space. In this section, 
an overview of a general-purpose, state-of-the-art numerical implementation is presented. 
Many of the features and techniques to be discussed, in this section, were developed previ- 
ously for elastostatics (e.g., Banerjee et al, 1985, 1988), and elastodynamics (e.g., Banerjee 
et al, 1986; Ahmad and Banerjee, 1988), but are here adapted for thermoelastic analysis. 

3.4.2 Temporal Discretization 

Consider, first, the time integrals represented in (3.2) as convolutions. Clearly, without 
any loss of precision, the time interval from zero to t can be divided into N equal increments 
of duration At. 

By assuming that the primary field variables, tp and up, are constant within each At 
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time increment, these quantities can be brought outside of the time integral. That is, 


N 

rnAt 

g 0a *b[X,t) = Y, t W 1 

1 90a(X-Z,t-T)dT 

n= 1 ^ 

(n— 1) At 

N 

rnAf 

90a * «/j(A, t) = ^ V/?(A) 

/ f0a(X-Z,t-r)dT 

l»=l 

/(n-l)At 


where the superscript on the generalized tractions and displacements, obviously, represents 
the time increment number. Notice, also, that, within an increment, these primary field 
variables are now functions of position only. Next, since the integrands remaining in 
(3.3) are known in explicit form from the fundamental solutions, the required temporal 
integration can be performed analytically, and written as 

/•nAt 

G N+i-n^ x -£)= g fia (X - e, t - r)dr (3.4a) 

P «/(n — 1) At 

y»nAt 

F * +1 -„ (x _ Q = / f pa [X - e,t - r)dr. (3 46) 

P J{n-l)At 

These kernel functions, and F£ a (X~t), are detailed in Appendix B.l. Combining 

(3.3) and (3.4) with (3.2) produces 

<*„(€)«7(0 = £ [ o*spo -<a +1 - n (^-o^wl^w- ( 3 - 5 ) 

n=l J * L J 

which is the boundary integral statement afte the application of the temporal discretiza- 
tion. 

3.4.3 Spatial Discretization 

With the use of generalized primary variables and the incorporation of a piecewise 
constant time stepping algorithm, the boundary integral equation (3.5) begins to show 
a strong resemblance to that of elastostatics, particularly for the initial time step (i.e., 
N = l). In this subsection, those similarities will be exploited to develop the spatial 
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discretization for the uncoupled quasistatic problem with two-dimensional geometry. This 
approximate spatial representation will, subsequently, permit numerical evaluation of the 
surface integrals appearing in (3.5). The techniques described here, actually, originated in 
the finite element literature, but were later applied to boundary elements by Lachat and 

Watson (1976). 

The process begins by subdividing the entire surface of the body into individual ele- 
ments of relatively simple shape. The geometry of each element is, then, completely defined 
by the coordinates of the nodal points and associated interpolation functions. That is, 

*(f)“*<(f) = tf-(f)*i- (36) 

with 

f intrinsic coordinates 
N w shape functions 
nodal coordinates 

and where u> is an integer varying from one to W, the number of geometric nodes in the 
element. Next, the same type of representation is used, within the element, to describe 
the primary variables. Thus, 



(3.7a) 


(3-76) 


in which u” u and t^ u are the nodal values of the generalized displacement and tractions, 
respectively, for time step n. Also, in (3.7), the integer « varies from one to fi, the total 
number of functional nodes in the element. From the above, note that the same number 
of nodes, and consequently shape functions, are not necessarily used to describe both the 
geometric and functional variations. Specifically, in the present work, the geometry is 
exclusively defined by quadratic shape functions. In two-dimensions, this requires the use 
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of three-noded line elements. On the other hand, the variation of the primary quantities 
can be described, within an element, by either quadratic or linear shape functions. (The 
introduction of linear variations proves computationally advantageous in some instances.) 

Once this spatial discretization has been accomplished and the body has been subdi- 
vided into M elements, the boundary integral equation can be rewritten as 

»».(«)<(?> = E { E [ [<„■"“"(*« - e)w„( f )i3„ 

n=l t m = 1 J l 

- F? a +l ~ n m) - e)^(o«?w]^(x(«:))} (3.8) 

where 

M 

S=^2 S m- 

m— 1 

In the above equation, t% u and v. n 0u are nodal quantities which can be brought outside the 
surface integrals. Thus, 

N ( M . 

«„uK (S) = E E - WMdsm;)) 

n = 1 t m= l Js m 

- f F»+ l - n {X[ S ) - ( 3 . 9 ) 

" Sm ' 

The positioning of the nodal primary variables outside the integrals is, of course, a key 
step since now the integrands contain only known functions. However, before discussing 
the techniques used to numerically evaluate these integrals, a brief discussion of the sin- 
gularities present in the kernels G^ a and Fg a is in order. 

The fundamental solutions to the uncoupled quasistatic problem contain singularities 
when the load point and field point coincide, that is, is when r = 0. The same is true of G” a 
and Fa , since these kernels are derived directly from the fundamental solutions. Series 

pCM. ' 

expansions of ter ms present in the evolution functions can be used to deduce the level of 
singularities existing in the kernels. 

A number of observations concerning the results of these expansions should be men- 
tioned. First, as would be expected F* p has a stronger level of singularity than does the 
corresponding G^p, since an additional derivative is involved in obtaining F*p from G a p. 
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Second, the coupling terms do not have as a high degree of singularity as do the corre- 
sponding non-coupling terms. Third, all of the kernel functions for the first time step could 
actually be rewritten as a sum of steady-state and transient components. That is, 

GU - * G a0 + tr Gi 0 
F l a0 =" F a0 + tr F* p . 

Then, the singularity is completely contained in the steady-state portion. Furthermore, 
the singularity in Gjy and F^ is precisely equal to that for elastostatics, while G eg and F gg 
singularities are identical to those for potential flow. (For two-dimensions, the subscript 
9 equals three.) This observation is critical in the numerical integration of the F a0 kernel 
to be discussed in the next subsection. However, from a physical standpoint, this means 
that, at any time t, the nearer one moves toward the load point, the closer the quasistatic 
response field corresponds with a steady-state field. Eventually, when the sampling and 
load points coincide, the quasistatic and steady-state responses are indistinguishable. As 
a final item, after careful examination of Appendix B.l, it is evident that the steady-state 
components in the kernels G" p and F£p, with n > 1, vanish. In that case, all that remains 
is a transient portion that contains no singularities. Thus, all singularities reside in the 
••G a0 and , ‘F a0 components of G l a& and F* p , respectively. 

3.4.4 Numerical Integration 

Having clarified the potential singularities present in the coupled kernels, it is now 
possible to consider the evaluation of the integrals in equation (3.9). That is, for any 
element m, the integrals 

f G N+1 - n (X({)-Z)N„U)dS(XU)) (3.10a) 

Js m 

f F N+1 ~ n (X(s) - £)N u (<)dS(X(s)) (3-106) 

Js m 

will be examined. To assist in this endeavor, the following three distinct categories can be 
identified. 
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(1) The point £ does not lie on the element m. 

(2) The point £ lies on the element m, but only non-singular or weakly singular integrals 

are involved. 

(3) The point £ lies on the element m, and the integral is strongly singular. 

In practical problems involving many elements, it is evident that most of the inte- 
gration occurring in equation (3.9) will be of the category (1) variety. In this case, the 
integrand is always non-singular, and standard Gaussian quadrature formulas can be em- 
ployed. Sophisticated error control routines are needed, however, to minimize the com- 
putational effort for a certain level of accuracy. This non-singular integration is the most 
expensive part of a boundary element analysis, and, consequently, must be optimized to 
achieve an efficient solution. In the present implementation, error estimates, based upon 
the work of Stroud and Secrest (1966), are employed to automatically select the proper 
order of the quadrature rule. Additionally, to improve accuracy in a cost-effective man- 
ner, a graded subdivision of the element is incorporated, especially when £ is nearby. For 
two-dimensional problems, the integration order varies from two to twelve, within each of 
up to four element subdivisions. 

Turning next to category (2), one finds that again Gaussian quadrature is applicable, 
however, a somewhat modified scheme must be utilized to evaluate the weakly singular 
integrals. This is accomplished in two-dimensional elements via suitable subsegmentation 
along the length of the element so that the product of shape function, Jacobian and kernel 
remains well behaved. 

Unfortunately, the remaining strongly singular integrals of category (3) exist only in 
the Cauchy principal value sense and cannot, in general, be evaluated numerically, with 
sufficient precision. It should be noted that this apparent stumbling block is limited to the 
strongly singular portions, -F ti and ••Fee, of the kernel. The remainder of F^, including 
tfpi an d tT Fgg, can be computed using the procedures outlined for category (2). However, 
as will be discussed in the next subsection, even category (3) t, F ij and ••Fee kernels can be 
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accurately determined by employing an indirect ‘rigid body’ method originally developed 
by Cruse (1974). 


3.4.5 Assembly 

The complete discretization of the boundary integral equation, in both time and space, 
has been described, along with the techniques required for numerical integration of the 


kernels. Now, a system of algebraic equations can be developed to permit the approximate 
solution of the original quasistatic problem. This is accomplished by systematically writing 
(3.9) at each global boundary node. The ensuing nodal collocation process, then, produces 
a global set of equations of the form 


N , 

£( 

n=l v 


>AT + l-n 


{> n- 


h^JV + l- 


1 - n ]{u B }) 


( 0 ), 


( 3 . 11 ) 


where 

j£.N+i-n] unassembled matrix of size ( d+ 1)P x (d+ 1 )Q, with coefficients determined 
from (3.10a) 

[jw+i-»] assembled matrix of size (d+ l)Px(<f+l)P, with coefficients determined from 
(3.10b) and c fia included in the diagonal blocks 

{t") global generalized nodal traction vector with (d+ 1 )Q components 
(u”} global generalized nodal displacement vector with (d + 1)P components 
{0} null vector with (d+ 1)P components 
P total number of global functional nodes 


Q — Em=l 

A m number of functional nodes in element m 
d dimensionality of the problem. 


In the above, recall that the terms generalized displacement and traction refer to the 
inclusion of the temperature and flux, respectively, as the (d+ 1) component at any point. 
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Consider, now, the first step. Thus, for N = 1, equation (3.11) becomes 

[G 1 ]{t 1 }-[J ?1 ]{u 1 } = {0>. (3-12) 

However, at this point the diagonal block of {F 1 } has not been completely determined 
due to the strongly singular nature of “Fa and •• F ee . Following Cruse (1974) and, later, 
Banerjee et al (1986) in elastodynamics, these diagonal contributions can be calculated 
indirectly by imposing a uniform ‘rigid body’ generalized displacement field on the same 
body, but under steady-state conditions. Then, obviously, the generalized tractions must 

be zero, and 

1"F]{1} = {0}, (3-13) 

where {1} is a vector symbolizing a unit uniform motion. Using (3.13), the desired diagonal 
blocks, **F i} and •'Fee, can be obtained from the summation of the off-diagonal terms of 
[** F). The remaining transient portion of the diagonal block is non-singular, and hence can 
be evaluated to any desired precision. With that step completed, (3.12) is rewritten as 

|C? 1 ]{t 1 }-[F 1 ]{u 1 } = {°>. (314) 

In a well-posed problem, at time At, the set of global generalized nodal displacements 
and tractions will contain exactly (d + l)P unknown components. Then, as the final stage 
in the assembly process, equation (3.14) can be rearranged to form 

( 315 > 


in which 

{s 1 } unknown components of {u 1 } and {t 1 } 
{y 1 } known components of {u 1 } and {t 1 } 
[A 1 ], [B 1 ] associated matrices 
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3.4.6 Solution 

To obtain a solution of (3.15) for the unknown nodal quantities, a decomposition 
of matrix [A 1 ] is required. In general, [ A 1 ] is a densely populated, unsymmetric matrix. 
The out-of-core solver, utilized here, was developed originally for elastostatics from the 
LINPACK software package (Dongarra et al, 1979) and operates on a submatrix level. 
Within each submatrix, Gaussian elimination with single pivoting reduces the block to 
upper triangular form. The final decomposed form of [A 1 ] is stored in a direct-access file 
for reuse in subsequent time steps. Backsubstitution then completes the determination of 
{z 1 }. Additional information on this solver is available in Banerjee et al (1985). 

After turning from the solver routines, the entire nodal response vectors, {u 1 } and 
{P}, at time At are known. For solutions at later times, a simple marching algorithm is 
employed. Thus, from (3.11) with N — 2, 

[G 1 ]^ 1 } - [F 1 ]!* 1 } + [G 1 ]^ 2 } - IF 1 ]!* 2 } = {0}. ( 316 ) 

Assuming that the same set of nodal components are unknown as in (3.14) for the first 
time step, equation (3.16) is reformulated as 

[A 1 ]!* 2 } = [iPKy 2 } - [G 2 ]^} + [F 2 ]!* 1 }. (3-17) 

Since, at this point, the right-hand side contains only known quantities, (3.17) can be 
solved for {z 2 }. However, the decomposed form of [A 1 ] already exists on a direct-access file, 
so only the relatively inexpensive backsubstitution phase is required for the solution. 

The generalization of (3.17) to any time step N is simply 

[A 1 ]!*"} = m{y N ) - E ([ G" +1 -"]{t" } - [F" +1 -"]{u"}) (3.18) 

n=l ' 

in which the summation represents the effect of past events. By systematically storing 
all of the matrices and nodal response vectors computed during the marching process, 
surprisingly little computing time is required at each new time step. In fact, for any time 
step beyond the first, the only major computational task is the integration needed to form 


16 



[G^J and [F^]. Even this process is somewhat simplified, since now the kernels are non- 
singular. Also, as time marches on, the effect of events that occurred during the first time 
step diminishes. Consequently, the terms containing [G*] and \F»] will eventually become 
insignificant compared to those associated with recent events. Once that point is reached, 
further integration is necessary, and a significant reduction in the computing effort per 
time step can be achieved. 

It should be emphasized that the entire boundary element method developed, in this 
section, has involved surface quantities exclusively. A complete solution to the well-posed 
linear uncoupled quasistatic problem, with homogeneous properties, can be obtained in 
terms of the nodal response vectors, without the need for any volume discretization. In 
many practical situations, however, additional information, such as, the temperature at 
interior locations or the stress at points on the boundary, is required. The next subsection 
discusses the calculations of these quantities. 

3.4.7 Interior Quantities 

Once equation (3.18) is solved, at any time step, the complete set of primary nodal 
quantities, {u n } and {t N }, is known. Subsequently, the response at points within the body 
can be calculated in a straightforward manner. For any point £ in the interior, the gener- 
alized displacement can be determined from (3.9) with cp a = 6p a . That is, 

<({) -£{£; k / c£. +I -"(*(c) - {)AU{)js(x( t )) 

n= 1 ^ m= 1 L 

~ j s F? a + 1 - n (X(<) ~ flhU?)dS(Xfc))] }. (3.19) 

Now, all the nodal variables on the right-hand side are known, and, as long as, £ is not on 
the boundary, the kernel functions in (3.19) remain non-singular. However, when £ is on the 
boundary, the strong singularity in **Fp a prohibits accurate evaluation of the generalized 
displacement via (3.19), and an alternate approach is required. The apparent dilemma is 
easily resolved by recalling that the variation of surface quantities is completely defined 
by the elemental shape functions. Thus, for boundary points, the desired relationship is 
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simply 


«?(0 = *«(?)“£ 


(3.20) 


where N w {<) are the shape functions for the appropriate element and ? are the intrinsic 
coordinates corresponding to ( within that element. Obviously, from (3.20), neither in- 
tegration nor the explicit contribution of past events are needed to evaluate generalized 

boundary displacements. 

In many problems, additional quantities, such a heat flux and stress, are also important. 
The boundary integral equation for heat flux, can be written 

N ( M r r 

«r(e) - £ { £ k / - e 

n=l t m= 1 JSm 


where 


EfciXis) -() = -k 


'0Oi 


e)iM?)ds(*(?))]}- 

(3.21) 

9G%,(X(s) - i) 

(3.21a) 

dii 

dF£ e (XU) - 0 

(3.216) 



This is valid for interior points, whereas, when £ is on the boundary, the shape functions 
can again be used. In this latter case, 


-Nw(?)<?w = n *(£)9, W (£) 


(3.22a) 




(3.226) 


which can be solved for boundary flux. Meanwhile, interior stresses can be evaluated from 


$(() = E { E k f - QN m (s)dS{Xto) 

n=l *- m=l L JSm 

- / o&'-Wf) - aw.b)js(jr(f))] } 

» 5 m 


(3.23) 


in which 


2 uv dG n 0t (dG%i 3G£- 

^ is) - f) = + “{sir + si 


- PSijG n $e 


(3.23a) 
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(3.236) 


2uv e 9F7 e fdF % dF%\ 

+lt {w + ~^r)~ 3 v 

Equation (3.23) is, of course, developed from (3.19). Since strong kernel singularities 
appear when (3.23) is written for boundary points, an alternate procedure is needed to 
determine surface stress. This alternate scheme exploits the interrelationships between 
generalized displacement, traction, and stress and is the straightforward extension of the 
technique typically used in elastostatic implementation (Cruse and Van Buren, 1971). 
Specifically, the following can be obtained 



»,■(*)«$(« = *«(?)*£ 

(3.24a) 

^ i}kl 

(<,(£) + «&(£)) = -0fciJMs)«2L 

(3.246) 


d S ^ as * w 

(3.24c) 


in which is obviously the nodal temperatures, and, 


D*j kl = A SijSki + 2 pSikSji- 

Equations (3.24) form an independent set that can be solved numerically for *£(« and 
u v.(£) completely in terms of known nodal quantities t& and tJL, without the need for 
kernel integration nor convolution. Notice, however, that shape function derivatives ap- 
pear in (3.24c), thus constraining the representation of stress on the surface element to 
something less than full quadratic variation. The interior stress kernel functions, defined 
by (3.23), are also detailed in Appendix B.l. 


3.4.8 Advanced Features 

The thermoelastic formulation has been implemented as a segment of the state-of-the- 
art, general purpose boundary element computer program, GP-BEST. Consequently, many 
additional features, beyond those detailed above, are available for the analysis of complex 
engineering problems. Perhaps, the most significant of these items, is the capability to 
analyze substructured problems. This, not only extends the analysis to bodies composed of 
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several different materials, but also often provides computational efficiencies. An individual 
substructure or geometric modeling region (GMR) must contain a single material. During 
the integration process, each GMR remains a separate entity. The GMR’s are then brought 
together at the assembly stage, where compatibility relationships are enforced on common 
boundaries between regions. Typically, compatibility ensures continuous displacement and 
temperature fields across an interface, however, recent enhancements to the code permit 
sliding between regions, spring contacts and interfacial thermal resistance to model air 
gaps or coating resistances. In the latter instances, discontinuities appear at the interface. 
In any case, the multi-GMR assembly process produces block-banded system matrices that 
are solved in an efficient manner. 

As another feature, a high degree of flexibility is provided for the specification of bound- 
ary conditions. In general, time-dependent values can be defined in either global or local 
coordinates. Not only can generalized displacements and tractions be specified, but also 
spring and convection boundary conditions are available. Another recent addition permits 
time-dependent ambient temperatures. A final item, worthy of note, is the availability of 
a comprehensive symmetry capability which includes provisions for both planar and cyclic 

symmetry. 

In the past year, an interface to the well-known PATRAN graphics package was de- 
veloped. This interface allows the user an option to view deformed shapes, temperatures 
and stress boundary profiles or contours. A number of PATRAN-produced illustrations 
are included in Section 5, however, in the next section, a couple of examples are presented 
to demonstrate the validity and applicability of this boundary-only formulation. 

3.5 Numerical Examples 

3.5.1 Sudden Heating of A luminum Block 

As a first example, transient heating of an aluminum block is examined under plane 
strain conditions. The block, shown in Figure 3.1, initially rests in thermodynamic equi- 
librium at zero temperature. Then, suddenly, the face at Y = 1.0 in. is elevated to 100°F, 
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while the remaining three faces are insulated and restrained against normal displacements. 
Thus, only axial deformation in the K-direction is permitted. Naturally, as the diffusive 
process progresses, temperature builds along with the lateral stresses »„ and «... To com- 
plete the specification of the problem, the following standard set of material properties are 

used to characterize the aluminum: 

E = 10 x 10 6 p«, v = 0.33, 

a = 13 x 10 ~ 6 /°F, 

k = 25in.lb./sec.in.°.F, pc, — 200in.lb./in. F. 

The two-dimensional boundary element idealization consists of the simple four element, 
eight node model included in Figure 3.1. A time step of 0.4 sec. is selected, corresponding 
to a non-dimensional time step of 0.5. Additionally, a finite element analysis of this same 
problem was conducted using a modified thermal version of the computer code CRISP 
(Gunn and Britto, 1984). The finite element model is also a two-dimensional plane strain 
representation, however, sixteen linear strain quadrilaterals are placed along the diffusion 
length. In the FE run, a time step of 0.2 sec. is employed. 

Temperatures, displacements, and stresses are compared in Table 3.1. Notice that the 
boundary element analysis, with only one element in the flow direction, produces a better 
time-temperature history than does a sixteen element FE analysis with a smaller time 
step. Both methods exhibit greatest error during the initial stages of the process. This is 
the result of the imposition of a sudden temperature change. Meanwhile, the comparison 
of the overall axial displacement indicates agreement to within 3% for the BE analysis 
and 5% for the FE run. A steady-state analysis via both methods produces the exact 
answer to three digit accuracy. The last comparison, in the table, involves lateral stresses 
at an integration point in the FE model. The boundary element results are quite good 
throughout the range, however, the FE stresses exhibit considerable error, particularly 
during the initial four seconds. Actually, these finite element stress variations are not 
unexpected in light of the errors present in the temperature and displacement response. 
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Recall that in the standard finite element process, stresses are computed on the basis of 
numerical differentiation of the displacements, whereas in boundary elements, the stresses 
at interior points are obtained directly from a discretized version of an exact integral 
equation. Consequently, the BE interior stress solution more nearly coincides with the 
actual response. 

3.5.2 Circular Disc 

Next, transient thermal stresses in a circular disc are investigated. The disc of radius 
‘o’ initially rests at zero uniform temperature. The top and bottom surfaces are thermally 
insulated, and all boundaries are completely free of mechanical constraint. Then, suddenly, 
at time zero, the temperature of the entire outer edge (i.e., r = o) is elevated to unity and, 
subsequently, maintained at that level. 

The boundary element model of the disc with unit radius is shown in Figure 3.2. Only 
four quadratic elements are employed, along with quarter symmetry. Ten interior points are 
also included strictly to monitor response. In addition, the following non-dimensionalized 
material properties are arbitrarily selected for the plane stress analysis: 

E = 1.333 pc, = 1.0 

v = 0.333 k = 1.0 

a = 0.75 

Results obtained under quasistatic conditions for a time step of 0.005 are compared, in 
Figures 3.3, 3.4 and 3.5, to the analytical solution presented in Timoshenko and Goodier 
(1970). Notice that temperatures, as well as radial and tangential stresses are accurately 
determined via the boundary element analysis. In particular from Figure 3.5, even the 
tangential stress on the outer edge is faithfully reproduced. 
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TABLE 3.1 

SUDDEN HEATING OF A CUBE 


Time 

(sec.) 

Temperature (°F) Axial Displacement (m in.) 

aty = o at y = 1.0 

Exact FE BEM Exact FE BEM 

Lateral Stress (ksi) 
at y = 0.5312 
Exact FE BEM 

0.8 

4.7 

3.4 

3.8 

910 

860 

920 

-5.6 -3.9 

-5.4 

1.6 

22.0 

19.8 

20.7 

1290 

1250 

1320 

-9.1 -7.7 

-9.2 

2.4 

38.3 

36.4 

37.7 

1570 

1540 

1610 

-11.3 -10.3 

-11.7 

3.2 

51.5 

50.0 

51.5 

1780 

1760 

1840 

-13.1 -12.2 

-13.5 

4.0 

61.9 

60.7 

62.2 

1950 

1930 

2000 

-14.4 -13.8 

-14.8 

4.8 

70.1 

69.1 

70.5 

2090 

2070 

2130 

-15.5 -15.0 

-15.9 

5.6 

76.5 

75.7 

76.9 

2200 

2180 

2230 

-16.3 -15.9 

-16.7 

6.4 

81.5 

80.9 

81.9 

2280 

2270 

2310 

-17.0 -16.7 

-17.3 

7.2 

85.5 

84.9 

85.8 

2340 

2330 

2370 

-17.5 -17.2 

-17.8 

8.0 

88.6 

88.2 

88.8 

2400 

2390 

2410 

-17.9 -17.7 

-18.1 
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4. INTEGRAL FORMULATION FOR FLUIDS 

4.1 Introduction 

Attention is now shifted to the hot fluid. A number of integral formulations will be 
presented for both incompressible and compressible thermoviscous flow. In particular, 
significant effort has been directed during the past year toward the development and im- 
plementation of the convective formulations. As a result, boundary element solutions can 
now be obtained in the high Reynolds number range. 

After a separation into the classes of incompressible and compressible flow, individual 
subsections present the governing equations, integral representations, numerical imple- 
mentation and numerical examples. It will be evident that a vast majority of the required 
work has been finished for the incompressible case. On the other hand, while compressible 
formulations are complete, most of the implementation effort is planned for 1990. 

4.2 Incompressible Thermoviscous Flow 

4.2.1 Introduction 

In the following, four distinct formulations are presented for incompressible flow: 
steady, time-dependent, steady convective, and time-dependent convective. The primary 
variables in each case are velocity, temperature, traction and heat flux. This is the set 
of variables for which boundary conditions are most readily defined, and for which the 
extension to three-dimensions is most easily accomplished. As will be seen, the individual 
formulations have much in common. The major differences involve the fundamental solu- 
tions that are employed, and the treatment of the contributions of past events. All four 
formulations are available within the computer code GP-BEST. 

4.2.2 Governing Equations 

Application of the Principles of the Conservation of Mass, Momentum and Energy for 
an incompressible thermoviscous fluid lead to the development of the following differential 
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equations: 


(4.1a) 


dvi 

dxi 

dp 


= 0 


9 2 t >i dp_ _ Dvi_ _ 

'dzjdxj dxi p Dt +I ' 

, d 2 e 


where 

Xj Eulerian coordinate 
t time 

^ velocity vector 
p pressure 
0 temperature 
p mass density 
n viscosity 

jfc thermal conductivity 
c f specific heat 
fi body force 
V> body source, 


(4.16) 

(4.1c) 


and the operator 


D _ d_ d 
Dt~ dt + Vj dxj 


(4.2) 


represents a material time derivative. By introducing a constant free stream velocity U, 
and a velocity perturbation such that 


v,- = Ui 4- Uj, 


(4.3) 


the governing equations can be rewritten as 


3ui 

d%i 


= 0 


d 2 t ii dp 

dxjdxj 


3ui TT dui duj , n 

- p-p- - pUj* P u fa t- fi - o 

dn dt 3 dxj dxj 


(4.4a) 


(4.46) 
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(4.4c) 


d 2 e 

' dxjdxj 


do TT do 30 

p^—. - P c * u >^r ~ pwzi : + ^ = °* 




dxj 


dxj 


Note that in equations (4.4) only the terms and pc t Uj-g~ are actually nonlinear, 

although in some instances the body forces and sources may also contain nonlinearities. A 
number of distinct integral formulations are possible, depending upon which of the linear 
terms are included in the differential operator. All terms excluded from the differential 


operator, axe then grouped together as effective body forces and sources, /' and t/>', respec- 
tively. Four particularly useful integral formulations are detailed in the next subsection. 


4.2.3 Integral Representations 
4.2.3. 1 Steady 

In this first formulation the time-dependent terms vanish, and the entire contribution 
of the convective terms are considered as effective body forces and sources. Thus, 


, dui dui 

X “ Mw, - '“’ S" + h 

, 98 96 

* = - pc ^'a^' pCeU ^ + v ’- 


(4.5a) 

(4.56) 


As a result, the well-known fundamental solutions for incompressible Stokes flow and 
steady-state heat conduction are applicable. The integral formulation, which can be de- 
rived directly from the governing differential equation (Dargush and Banerjee, 1990b), can 


be written 


CapUa ~ [ [G a pt a — F a pU a — G a pt°]dS + [ [Dapk&ka + G a fifa ] dV 
J S Jv 


(4.6) 


where 


U a = («1 “2 6} ( 4,7 °) 

t a = {ti *2 q} ( 4 - 7fc ) 

fa ~ {/l h ( 4,7c ) 

are generalized velocities, tractions, and body forces. In (4.7b), U are the surface tractions 

defined by 

ti = r.yny - pn, (4.8a) 
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with rii representing the local unit outward normal to the surface 5, and r,y the fluid 
stresses, while the heat flux is defined via 



(4.86) 


Furthermore, 


— 


Cij 0 
0 C0o 


G a p — 


Ga 


0 

0 Ggff 


F a p = 


Fa 0 
0 Fg$ 


D n Rk = 


dG a p 


' afik dx k 

a ka = \p{Gk + PCM + »k)0] 


t° a = °ka n k- 


(4.9a, 6, c) 


(4.9<f) 

(4.10a) 

(4.106) 


In the terminology of Lighthill (1952), is the momentum flux tensor or fluctuating 
Reynolds stress. Here, a£ a is labeled the generalized convective stress tensor, while t° a is 
the generalized convective traction. Both o° ka and t° a contain terms which are nonlinear in 
the generalized velocities. 

In (4.9a), c, y (£) and <**(£) are constants. When i is inside S,c i} = and c ee = l. If 
£ is on the boundary then the values are determined by the relative smoothness of 5 at 
£. For £ outside the region V, both c i} and cg$ are zero. Meanwhile, the kernel functions 
Gij,G ee ,Fij and F ee are provided in Appendix B.2. 


4.2.3. 2 Time-Dependent 

For this next formulation, the effective body forces and sources are identical to those 
provided in (4.5), however, the time-dependent terms are now included in the linear oper- 
ator. The required fundamental solution for the viscous portion was first given by Oseen 
(1927), while the transient heat conduction fundamental solution is well-known (Carslaw 
and Jaeger, 1959). By applying standard methodology (Banerjee and Butterfield, 1981; 
Dargush and Banerjee, 1990c), the following governing integral equations can be derived 

CapVia = [ [gap * - fap * «« ~ 9aP * C] dS + f [d a p k * °ka + 9a fi * fa ~ 9afiP<] & ( 4 ' U ) 

is Jv 
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Note that (4.11) is similar to (4.6) for the steady case, except that Riemann convolution 
integrals over time have been introduced, along with an initial condition volume integral 
involving u%. Kernel functions, G a0 and F a0 , developed from the instantaneous point force 
and source adjoint fundamental solutions g a0 and are provided in Appendix B.3. 


4.2.3. 3 Steady Convective 

At large Reynolds number the integral formulations presented in the previous two 
subsections are not suitable, because the nonlinear convective terms involving t° a and o° ka 
dominate the problem. For high speed flow, it is essential to include more of the physics 
of the problem in the fundamental solution. This can be accomplished by including the 
linear convective terms P U j and pc ( U,£- in the differential operator. Then the effective 


body forces and sources become 


/; - + u 


(4.12a) 


and 

i/>' = -pc t Uj~- + r/) (4.126) 

cf Xy 

respectively. The corresponding integral equations, under steady conditions, are simply 


«o«. = / [O. - f.V. - <» + jf dv - < 413 > 

The superscript U on the kernel functions is a reminder that these are based upon convec- 
tive fundamental solutions. The kernel G u a0 is detailed in Appendix B.4. Meanwhile, 


o k ° = [pufctti pc t u k 6\ 


(4.14a) 


*U° — rr U °„. 

*a =Vka n k- 


(4.146) 


Interestingly, both the convective fundamental solution and integral equation for viscous 
flow were developed by Oseen in the early portion of this century. In fact, formulations 
similar to (4.6), (4.11) and (4.13) are all presented in an elegant manner by Oseen in his 
1927 monograph. Of course, this was well before the advent of the computer. As a result, 
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Oseen was unable to do much with his formulations other than some approximations at 
very low Reynolds number. 

At first glance, equations (4.6) and (4.13) are quite similar. However, the significant 
differences in the behavior of the kernel functions necessitate quite different numerical 
treatment. The exact nature of these kernels will be examined later during the discussion 
of numerical integration. 

4.2.3.4 Time-Dependent Convective 

The integral equations for this final case can be written formally as 

c a0 u a = f [9a0 * “ fafi * U <* ~ 9a0 * *a°] 

Js 

+ l [d l U * *£? + C * /- - dV ( 415 ) 

However, the instantaneous point force and source adjoint fundamental solutions g ^ cannot 
be easily expressed in terms of recognized mathematical functions. As an alternative, 
efforts are underway to develop polynomial approximations for g^ p over selected ranges of 
the independent parameters. 

4.2.4 Numerical Implementation 

4.2.4. 1 Introduction 

Analytical solutions are possible for only the simplest geometries and boundary con- 
ditions. More generally, approximations must be introduced in both time and space to 
expose the practical utility of these integral equations. Consequently, in this section, state- 
of-the-art boundary element technology is applied to steady and unsteady incompressible 
thermoviscous flows. Recent boundary element developments in the fields of elastodynam- 
ics (Banerjee et al, 1986; Ahmad and Banerjee, 1988) and thermoelasticity (Dargush and 
Banerjee, 1989b, 1990a) are directly applicable for these problems. The presentation below 
will concentrate on those aspects of the numerical implementation which differ from that 
detailed in Section 3. The current implementation is limited to the two-dimensional case, 
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although certainly all of the integral formulations presented in the previous subsection are 
equally valid in three dimension. 


4.2.4.2 Temporal and Spatial Discretization 

For time-dependent problems, the total time interval from zero to r is subdivided into 
N equal increments of duration At. Then, the field variables t a ,u a ,t£, and cr£ a are assumed 
constant within each Ar time increment. As a result, 


N rnAr ^ 

g a0 . t. « E « / 9 °» dt = £ K<%r +X 

n=l „=1 


( 4 . 16 ) 


with similar expressions holding for the remaining convolution integrals. This is identical 
to the treatment discussed in Section 3 for thermoelasticity. 

The methodology employed for spatial discretization of the bounding surface also fol- 
lows that described in Section 3. Thus, quadratic or linear shape functions are utilized to 
portray the functional behavior of the field variables over three-noded surface elements. 

However, in addition to the surface description, the domain must be discretized into 
cells in the regions where the nonlinear convective effects are important, or where nonzero 
initial conditions are present. Shape functions are once again introduced to approximate 
the geometric and functional variation with each volume cell. Thus, for any point X within 
an individual cell 

x<(?) = M w ($)x iw ( 4 - 17 ) 


*£,(?) = ( 418 ) 

where 

M W ,M U shape functions 
x iw nodal coordinates 
a? nodal generalized convective stress . 
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The current implementation utilizes six and eight-noded cells for the geometric represen- 
tation, along with both linear or quadratic functional variation. Typical cells are depicted 

in Figure 4.1. 

As a result of the spatial discretization, the boundary integral equation for time- 


dependent thermoviscous flow can now be written 


= £ I E [c / c»r+'N„is - / s _ - c / s _ G .V" + ‘ W "‘ ,S 

n= 1 m= 1 m 

+ 1 [«. / v o^- +i } + t p.. / Vi ■«'] 


(4.19a) 


while for steady conditions this reduces to 

c a 0 u a = 5Z [*<*« f G a pN u dS — U au f F a pN w dS — t au f G a pN^ 
m^l L «« • /s "> Sm 

+ E f D a 0kM u dV 

Jv ‘ 


dS 


(4.196) 


where M and L are the total number of surface elements and volume cells, respectively, 
and 

5 = £ (4.20a) 

m= 1 

V = ^Vj. ( 4 - 20fc ) 

/=1 

The positioning of the nodal variables outside of the integrals is a key step, since now the 
integrands of (4.19) contain only known functions, which can be evaluated numerically. 

Up to this juncture, the region of interest has been assumed to be composed of a single 
volume V with surface S. However, this need not be the case. In general, space may 
be subdivided into a number of individual non-overlapping geometric modeling regions 
(GMRs). Each GMR occupies a certain volume of space, say V 0 , bounded by the surface 
S g . For a point £ within V g , the integration required by (4.19) need only be conducted over 
S g and V g , since the contribution to u a (0 from the other GMRs outside S g will be zero. 
As a result, integration costs can be dramatically reduced by introducing multiple GMRs 
for thermoviscous flow problems. Additionally, there is no inherent requirement that all 
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GMRs utilize the same physical model. For example, one GMR could employ the steady 
formulation of equation (4.6), while a second region includes the convective kernel effects 
contained in the formulation of (4.13). In any case, compatibility must, of course, be 
maintained across all GMR-to-GMR interfaces. Examples of mixed GMR formulation are 
contained in Section 4.2.5 and form the basis of the approach for fluid structure interaction 
that will be explored in Section 5. 

4. 2.4. 3 Integration 

The evaluation of the integrals appearing in (4.19) is the next process to be examined. 
Due to the singular nature of the kernel functions G a p, F a p and D a pk considerable care must 
be exercised during numerical integration. This is particularly true for incompressible 
viscous flow, in which the final solution is extremely sensitive to errors in integration 
coefficients. In general, the integration algorithms must be much more sophisticated than 
those developed for thermoelasticity. In the present implementation, discussed in detail 
in Honkala and Dargush (1990), a number of different integration schemes are employed 
depending upon the order of the kernel singularity, the proximity of the field point £ to 
the element, and the size of the element. 

Before discussing the techniques used for numerical integration it is instructive to 
examine the nature of the kernel functions that are to be integrated. In particular, the 
convective kernels are quite different from those of elasticity and potential flow upon which 
most of the boundary element integration algorithms have been based. Consider first the 
two-dimensional infinite space response to a point force in an incompressible viscous fluid 
moving with a uniform reference velocity U\. The response Gn due to a unit force at (£i, £ 2 ) 
in the indirection is displayed in Figure 4.2a for points along i 2 = £2 at several values of 
Ui. Notice that as U\ increases the response becomes much more localized. At very high 
speed, Gn varies sharply only in a small band about the load point. Also, with non-zero 
Ui, the response is no longer symmetric in the indirection. This kernel contains a Doppler 
effect. The thermal portion of the kernel, G 33 , has a similar behavior and is shown in Figure 
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4.2c. However, Figure 4.2b depicts G 22 which is even more complicated, containing a sign 
reversal and inflection point near the load. (The remaining terms, G 12 ,G 13 ,G 2 i and G 23 
are zero along x 2 = £ 2 .) The traditional integration algorithms are not able to accurately 
capture this extremely localized behavior at large Reynolds number, and as a result new 
schemes have been devised to permit solutions for high speed flows. Details are provided 
below. 

Once again consider the following three distinct categories for the surface integrals: 

(1) The point £ does not lie on the element m. 

(2) The point £ lies on the element m, but the kernels involve only weakly singular inte- 

grands of the In r type. 

(3) The point £ lies on the element m, and the integral has a strong £ singularity. 

In practical problems involving many elements, it is evident that most of the integration 
occurring in equation (4.19) will be of the Category (l) variety. The integrand is non- 
singular and standard Gaussian quadrature can be employed. However, for near-singular 
cases when £ is close to element m very high order formulas are needed to capture the kernel 
behavior. For these instances, it is beneficial to identify the point X° on the element nearest 
to £, and then subdivide the interval of integration about X°. For non-convective kernels, 
within each of the two subsegments a nonlinear transformation is used to further reduce 
the order of Gaussian quadrature needed for high precision. This nonlinear transformation 
is similar to that proposed by Mustoe (1984) and Telles (1987), however it should be 
emphasized that subsegmentation is still required. For the convective near-singular case, 
graded subsegmentation is employed about X°. The smallest subsegments utilized are 
.0001 times the element length. High order Gauss formulas are used in segments near X°, 
while lower order formulas are used elsewhere. 

Turning next to Category (2), one finds that, unlike elasticity or potential flow, stan- 
dard Gaussian formulas alone are inadequate. Instead the terms involving In r must be 
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isolated and integrated with special log-weighted Gaussian integration. The remaining non- 
singular terms comprising G aP are then evaluated utilizing standard quadrature. Heavy 
subsegmentation is again included for convective kernels. 

The strongly singular integrals of Category (3) exist only in the Cauchy principal 
value sense and cannot be evaluated numerically with sufficient precision. Fortunately, 
the indirect ‘rigid body’ or ‘equipotential’ method, originally developed by Cruse (1974), 
is applicable, and leads to the accurate determination of the singular block of the second 
integral in (4.19). The remainder of that integral is non-singular. Consequently, subseg- 
mentation along with standard Gaussian quadrature is adequate. 

Similar care is needed for the volume integrals, which involve the kernel D afik con- 
taining a --type singularity. However, for two-dimensional volume integration, this kernel 
is only weakly singular, and can be evaluated in the following direct manner. First, the 
nearest node, say A, in cell l to the point $ is determined. The cell is then subdivided into 
triangles radiating from A a s shown in Figure 4.3. Next, each triangle is mapped onto a 
unit square. The apex corresponding to A is stretched to form one side of the square. This 
process essentially eliminates the 1 singularity. Finally, the square is further subsegmented 
in both radial and tangential directions depending upon the closeness of £ and the size of 
cell /. Standard Gaussian quadrature is applied to each subsegment. This cell integration 
scheme was based on work by Mustoe (1984) for elastoplasticity. In the present incom- 
pressible viscous flow implementation, tolerances have been tightened so that additional 
subsegmentation is performed, along with higher order quadrature formulas. For convec- 
tive kernels, the subsegmentation required is much more intense, and much higher order 
Gauss formulas are employed in the vicinity of the singularity. 

In time-dependent problems, beyond the first time step, additional integration is re- 
quired. This integration involves the kernels G^ p ,F2 P and D^ pk for n > 1. From Table 4.1, 
these are all nonsingular. As a result, a much less sophisticated integration scheme is em- 
ployed to obtain the required level of accuracy with fewer subsegments and gauss points. 
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If the initial velocities are not uniform, then the nonsingular initial condition integral of 
equation (4.19a) must also be evaluated at each time step. 


Table 4.1 - Kernel Singularities 


Kernel 
G n a0 for n > 1 

*u 

Fa 0 for n > 1 

D l a0k 

Kpk for n > 1 


Singularity Order 
In r 

non-singular 

i 

r 

non-singular 

i 

r 

non-singular 


4.2.4.4 Assembly 

Once the spatial discretization and numerical integration algorithms are completely 
defined, a system of nonlinear algebraic equations can be developed to permit an approx- 
imate solution of the thermoviscous boundary value problem. The method of collocation 

is employed by writing (4.19) at each functional mode. 

For each time step AT of a transient problem, this nodal collocation process yields 

y-' J G N-n+l t n _ F AT-n+l u n _ qN-u+ + j _ r ^ u « _ 0 (4.21) 

n= 1 

where 

t n nodal traction vector for time step n with 3Q components 

u « nodal velocity vector for time step n with 3P components 

t on Nodal convective traction vector for time step n with 3Q components 

& on nodal convective stress vector for time step n with 6P components 

u o nodal initial velocity vector with 3P components 
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G n 


D n 


T" 


M 

Q= £ A« 

m= 1 


unassembled matrix of size 3P x 3Q calculated from the first 
integral of (4.19) during time step n 

assembled matrix of size 3P x 3P calculated from the second 
integral of (4.19) during time step n, plus the c a p contribution 
in F 1 

assembled matrix of size 3P x 6P calculated from the first volume 
integral of (4.19) 

assembled matrix of size 3P x 3P calculated from the initial condition 

integral of (4.19) 

total number of functional nodes 

number of functional nodes in element m . 


All of the coefficient matrices in (4.21) contain independent blocks for each GMR in mul- 
tiregion problems. However, for any well-posed problem, the boundary conditions and 
interface relations remove all but 3P unknown components of u" and t". Furthermore, 
by solving (4.21) at each increment of time, all of the components of u n ,t n ,t on and c on for 
n < N are known from previous time steps. Then, (4.21) can be rewritten at time NAt as 

g(x) = Ax" - D 1 ^ + GH°" - By" 

N-l 

- X) [G" -n+1 t n - F" -n+1 u’ 1 - G" -n+1 t on + D" -n+1 <7 <m ] + r"u° = 0 (4.22) 


n= 1 


in which 




..N 


nodal vector of unknowns with 3P components 
nodal vector of knowns with 3Q components 


while A and B are the associated coefficient obtained from F 1 and G 1 . The A matrix 
now includes the compatibility relationships enforced on GMR interfaces. As a result, the 
GMR blocks in A are no longer independent, however A does remain block banded. 
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The terms included in the summation of (4.22) represent the contribution of past 
events. This, along with the terms By N and T N u°, can be simply evaluated once at each 
time step N with no need for iteration. Let, 

b N _ _ By N _ jf* 1 | G W-n+l t n _ jW-n+1^ _ QV-r,+ l t on + jjW-n+l^onj + r W u°. (4.23) 

n= 1 

Then (4.22) becomes the following nonlinear set of algebraic equations 

g(x) = Ax n - D 1 a oN + G 1 t oJV ’ + b N = 0. (4-24) 


A closer examination of b A is in order. For example with N — 1 


b 1 = -By 1 +r 1 u°, 


(4.25a) 


while for the second time step 

b 2 = -By 2 - G 2 t 1 + F 2 ^ + G 2 t o1 - D 2 ff o1 + T 2 u° (4.256) 

Obviously, for each step N, one new set of matrices G^.F^.D" and T N must be determined 
via integration and assembly. Integration, particularly the volume integration needed for 
D n and T N , can be quite expensive. 

As an alternative to the convolution approach defined above, a time marching recur- 
ring initial condition algorithm can be employed. This has been utilized by a number of 
researchers for transient problems of heat conduction, acoustics, and elasticity (Banerjee 
and Butterfield, 1981). For this latter approach, at time step N the entire contribution of 
past events is represented by an initial condition integral which utilizes u"" 1 as the initial 

velocity. Thus, 

g(x) = Ax n - DV n + G x t oW + b N = 0 (4-26) 

with 

b N = -By N + rV*- 1 . ( 4 - 27 ) 
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Obviously, (4.26) is identical to (4.24). Only the evaluation of b" is different. The advan- 
tage of the recurring initial condition approach is that no integration is needed beyond the 
first time step. However, volume integration is required throughout the entire domain be- 
cause of the presence of u"" 1 , even for linear problems in which volume integration would 

not normally be required. 

In order to take full advantage of both methods, the present work utilizes the con- 
volution approach in linear regions, and the recurring initial condition algorithm for the 
remaining nonlinear GMRs which are filled with volume cells. Since b N can be computed 
independently for each GMR, this new dual approach provides no particular difficulty. 


4.2.4. 5 Solution 


An iterative algorithm, along the lines of those traditionally used for BEM elastoplas- 
ticity (Banerjee and Butterfield, 1981; Banerjee et al, 1987), can be employed to solve the 
boundary value problem. However, convergence is usually achieved only at low Reynolds 
number. More generally the interior equations must be brought into the system matrix, as 


in (4.24), and a full or modified Newton- Raphson algorithm must be employed to obtain 
solutions even at moderate Reynolds number. (Similar 'variable stiffness algorithms have 
also been introduced by Banerjee and Raveendra (1987) and Henry and Banerjee (1988) 
for elastoplasticity.) Symbolically, at any iteration fc, 


dg 

dx 


|Ax fc | = -|g(x) fc | 


( 4 . 28 ) 


where 


x fc+1 = x fc + Ax fc 


( 4 . 29 ) 


and the derivatives on the lefthand side of (4.28) axe evaluated at x fc . With the full Newton- 
Raphson approach, the system matrix must be formed and decomposed at each iteration. 
The out-of-core solver used in the present implementation was developed originally for 
elastostatics (Banerjee et al, 1985) from the LINPACK software package (Dongarra et al, 
1979), and operates on a submatrix level. Within each submatrix, Gaussian elimination 
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with single pivoting reduces the block to upper triangular form. The final decomposed 
compacted form of the system matrix is stored in a direct access file for later reuse. Back- 
substitution completes the determination of Ax fc . Iteration continues until 


H(Ax w ni 

l!(x") k |l 


(4.30) 


where « is a small tolerance, and ||x|| is the Euclidean norm of x. For the modified Newton- 
Raphson algorithm, the system matrix is not formed at every iteration, and only backsub- 
stitution is needed to determine Ax k . 


4.2.4.6 Calculation of Additional Boundary Qua ntities 

Once the iterative process has converged, a number of additional boundary quantities 
of interest can be easily calculated. For example, lift and drag can be calculated by numer- 
ically integrating the known nodal traction and shape function products over the surface 
elements of interest. Low order Gaussian quadrature is adequate for this integration, since 
all the functions are very well behaved. 

Furthermore, at each boundary node, the pressure p, stress <r,y, and strain rates ^ 
can be determined by simultaneously solving the following relationships: 




(4.31a) 


^ ~ M (S (0 + S (0 ) + P( ° " ° 

dxj dvn . . dN u 


(4.31fc) 

(4.31c) 

(31d) 


It should be emphasized that (4.31) represents a set of nine independent equations which 
are written at the boundary point f, and can be solved easily for p,<r i} and at that 
point. Afterward, boundary vorticity and dilatation can be obtained, respectively, from 


3u2 3ui 
dxi dx2 


(4.32a) 
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( 4 . 326 ) 


. du\ du2 

A = — - H 

dxi dx 2 

Of course, for incompressible flow, the dilatation should be zero, but (4.32b) can be used 
as a check. 

A comprehensive PATRAN interface has also been developed. Consequently, any of 
the quantities computed above may be displayed graphically in the form of profiles or 
contours. 

4.2.5 Numerical Examples 

4. 2. 5.1 Introduction 

All of the formulations discussed above have been implemented as a segment of GP- 
BEST, a general purpose boundary element code. In this section, a number of examples 
are included, primarily, to demonstrate the validity and attractiveness of the boundary 
element formulations. 

4. 2. 5.2 Converging Channel 

The two-dimensional incompressible flow through a converging channel also possesses 
a well known analytical solution which is purely radial (Millsaps and Pohlhausen, 1953). 
A comprehensive finite element study of this problem has been made by Gartling et al 
(1977). 

The boundary element model is shown in Figure 4.4a. The mesh contains 96 cells 
and is divided into two regions. The boundary conditions were modeled using an exact 
specification of the boundary conditions appearing in the analytical solution (Fig. 4.4a). 
Viscosity is unity, and tractions and density are incremented to reach higher Reynolds 
numbers. The Reynolds number for this problem is defined as 

(4.33) 

where F 2 (Rj) is the maximum velocity in the region, which is -24.0 for the problem solved 
here. 
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Figure 4.4b illustrates the results for two Reynolds numbers, indicating good accuracy 
along the entire width of the channel. Not only are the velocities accurate, but the pressures 

and tractions are very accurate also. 

It has been observed that finite element versions of this problem have several pecu- 
liarities which prevent the analytical solution from being reproduced. First of all, since 
velocities are often specified at the inlet and at the wall and centerline, ambiguous bound- 
ary condition specification results. Also, typically a parabolic “fully developed” velocity 
profile is usually specified at the inlet. However, the nonlinear solution has a flattened 
velocity distribution across the width of the channel (see Fig. 4.4b). Hence, the analyt- 
ical solution cannot be reproduced exactly if the “fully developed” profile is specified at 
the inlet. Also, the finite element modelers of this problem usually leave out the traction 
distribution at the exit and specify zero tractions there. This also gives rise to non-radial 

flow. 

The reason for so much interest in the converging flow problem is that it is one of 
the few problems possessing an analytical solution. However, by specifying a model which 
does not correspond to this problem, as in the finite element case, one cannot accurately 
compare results to the analytical solution. Any such comparisons are merely qualitative. 
In this light, the boundary element model here has utilized an exact model of the boundary 
condition and a meaningful comparison can be made. 

4.2.5. 3 Transient Couette Flow 

Consider as the first transient analysis the case of developing Couette flow between 
two plates, parallel to the x-z plane, a distance h apart. Initially, both of the plates, as 
well as the fluid, are at rest. Then, beginning at time t = 0, the bottom plate is moved 
continuously with velocity V in the x-direction. Due to the no-slip condition at the fluid- 
plate interface, Couette flow begins to develop as the vorticity diffuses. Eventually, when 
steady conditions prevail, the x-component of the velocity assumes a linear profile. 
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The following exact solution to this unsteady problem is provided by Schlicting (1955): 
v x {y, t) = V | er f c \ 2nT )i - 

ln=0 


where 


V = 


(4 a it/py / 2 


- Y^erfc[2(n+l)in - rj] 1 

(4.34a) 

n= 0 ) 


= 0 

(4.346) 

h 

ni (4a it/p) 1 ' 2 

(4.35a, 6) 


(4.35c) 

, and dv x /dx are zero. 



The two-dimensional boundary element model, utilized for this problem, is displayed 
in Figure 4.5. Four quadratic surface elements are employed, with one along each edge 
of the domain. A number of sampling points are included strictly to monitor response. 
Notice that the region of interest is arbitrarily truncated at the planes x = 0 and x = t. All 
of the boundary conditions are also shown in Figure 4.5. For the presentation of GPBEST 
results, all quantities are normalized. Thus, 


Y = ^ (4.36a) 

n 

T=~ ( 4 . 366 ) 

and the horizontal velocity is v x /V. Figure 4.6 provides the velocity profiles at four different 
times, using a time step AT = 0.025 and the convolution approach. There is some error 
present at small times near the top plate, where the velocity is nearly zero. Results at 
Y = 0.5 versus time are shown in Figure 4.7 for several values of the time step. Obviously, 
the correlation improves with a reduction in time step and AT = 0.025 provides accurate 
velocities throughout the time history. However, even for a very large time step, the 
GPBEST solution shows no signs of instability. Error, evident in the initial portion, 
diminishes with time, and all values of AT produce the correct steady response. Further 
reduction of AT beyond 0.025 yields little benefit. Instead, mesh refinement in the y- 
direction is needed, primarily to capture the short time behavior. Figure 4.8 shows the 
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GPBEST results for a model with just two, equal length, elements along each vertical side. 
The correlation with the analytical solution is now excellent. The time step selected for 
the refined model was based upon the general recommendation that 

AT a 0 05£ ”»" , (4.37) 

where is the length of the smallest element. 

The convolution approach, defined by equation (4.22), was used to obtain the results 
presented in Figures 4.6-4. 8. Alternatively, the recurring initial condition algorithm can 
be invoked. In that case, complete volume discretization is required even for this linear 
problem. For the model of Figure 4.6, a single volume cell connecting the eight nodes is 
all that is required. The GPBEST results for different values of AT are shown in Figure 
4.9. The solutions are good for the two smaller time step magnitudes, however there is a 
slight degradation in accuracy from the convolution results. 

Interestingly, the solution in (4.34a) is identical to that for one-dimensional transient 
heat conduction in an insulated rod with one end maintained at temperature V , while the 
other remains at zero. However, in a corresponding boundary element analysis, the numer- 
ical integrations defined in (4.19a) must be calculated much more precisely for unsteady 
viscous flow than for heat conduction in order to obtain comparable levels of accuracy. 


4.2.5.4 Flow Between Rotating Cylinders 

As the next example, the developing flow between rotating cylinders is analyzed. The 
inner cylinder of radius r< is stationary, while the outer concentric cylinder with radius 
r 0 is given a tangential velocity V, beginning abruptly at time zero. The steady solution 
appears in Schlicting (1955). However, even for the transient case, the flow is purely 
circumferential. Thus, the governing Navier-Stokes equations reduce to 


1 dvg V 0 \ dv e _ 

r dr r 2 ) ^ dt 

(4.38a) 

1 

Qj ! Qj 
1^3 

+ 

^ |<?M 

II 

o 

(4.386) 


46 



in polar coordinates ( r,6,z ). As discussed in Batchelor (1967), separation of variables can 
be used to obtain the following solution (Honkala and Dargush, 1990) 


where 


ve 


v r [r,t) = 0 

(r, t) = c 1 r+ C ^+f^D n {J 1 (X n r)Y 1 {X n r 0 )-Y 1 (X n r)J 1 {X n r 0 )} t 


,-K<* 


n= 1 


Cl = 


Fr c 
ri — r? 


C 3 = — Cifj 


Fin = — Ci[r„ J2{Xn r o) ~ fi^MAn**)] + C2[Jo(^n»V>) — *fo(A n f\)] 

F 2n = 0 i\rlY 2 (X n r o ) — r^Y 2 (X n ri)] — C 2 [Y 0 (X n r 0 ) — yi)(^n r *)] 


and A n is the nth root of the equation 


(4.39a) 

(4.396) 


(4.40a, b) 
(4.40c) 
(4.40(f) 
(4.40e) 


^(Ar^Y^Aro) - Ji(Ar 0 )yi(Ari) = 0. 


(4.41) 


Figure 4.10 depicts the boundary element model representing the region between the 
two cylinders. A thirty degree segment is isolated, with cyclic symmetry boundary condi- 
tions imposed along the edges 6 = 0° and 9 = 30°. The inner radius is unity, while an outer 
radius of two is assumed. Unit values are also taken for the viscosity, density and V . The 
model consists of six quadratic elements and two quadratic cells. The cells, of course, are 
not needed for linear analysis utilizing the convolution approach. 

Results of the GPBEST analysis are compared to the exact solution in Figure 4.11 
for convolution and in Figure 4.12 for the recurring initial condition algorithm. In both 
diagrams, results with and without the nonlinear convective terms are plotted. The re- 
sults are quite good throughout the time history with the convolution approach, while 
some noticeable error is present at early times for the recurring initial condition solutions. 
The linear and nonlinear velocity profiles are nearly identical, as expected from the exact 
solution expressed in (4.39b). However, unlike the previous example, the nonlinear terms 
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do not simply vanish from the integral equation written in cartesian form. Instead, the 
nonlinear surface and volume integrals must combine in the proper manner to produce 
the correct solution. Consequently, this problem provides a good test for the entire BEM 
formulation. 

Relative run times are shown in Table 4.2 for the different analysis types. Obviously, 
the nonlinear convolution approach is very expensive, since this involves volume integration 
at each time step. As a result, in the general implementation, convolution is only utilized 
in linear GMRs. 


Table 4.2 - Flow Between Rotating Cylinders 
(Run Time Comparisons) 


Analysis Type Time Marching Algorithm Relative CPU Time 


Linear Convolution 1.0 

Nonlinear Convolution 25.8 

Linear Recurring Initial Condition 1.5 

Nonlinear Recurring Initial Condition 1.8 


4.2.5. 5 Driven Cavity Flow 

The two-dimensional driven cavity has become the standard test problem for incom- 
pressible computational fluid dynamics codes. In a way, this is unfortunate because of the 
ambiguities in the specification of the boundary conditions. However, numerous results 
are available for comparison purposes. 

The incompressible fluid of uniform viscosity is confined within a unit square region. 
The fluid velocities on the left, right and bottom sides are fixed at zero, while a uniform 
nonzero velocity is specified in the x-direction along the top edge. Thus, in the top corners, 
the x-velocity is not clearly defined. To alleviate this difficulty in the present analysis, the 
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magnitude of this velocity component is tapered to zero at the corners. 

Results are presented for the four region, 324 cell boundary element model shown in 
Figure 4.13. Notice that a higher level of refinement is used near the edges. Spatial plots 
of the resulting velocity vectors are displayed in Figures 4.14a and b for Reynolds numbers 
(Re) of 400 and 1000, respectively. Notice that, in particular, the shift of the vortical 
center follows that described by Burggraf (1966) in his classic paper. A more quantitative 
examination of the results can be found in Figure 4.15 where the horizontal velocities on 
the vertical centerline obtained from the present GPBEST analysis are compared to those 
of Ghia et al (1982). It is assumed that the latter solutions are quite accurate since the 
authors employed a 129 by 129 finite difference grid. As is apparent, from the figure, all 
of the solutions are in excellent agreement. Finally, it should be noted that the simple 
iterative algorithm fails to converge much beyond Re = 100. Beyond that range the use of 
a Newton-Raphson type algorithm is imperative. 

4.2.5.6 Transient Driven Cavity F low 

The next example involves the initiation of flow in the same square cavity. An in- 
compressible fluid of uniform density and viscosity is at rest within a unit square region. 
The velocities of the vertical sides and the bottom are fixed at zero throughout time. At 
time zero, the horizontal velocity of the top edge is suddenly raised to a value of 1000 
and maintained at that level. A gradual transition of velocities is introduced near the top 

corners to provide continuity. 

The four region, 324 cell model shown in Figure 4.13 is employed for the boundary 
element analysis. The resulting velocity vector plots at several times are shown m Figure 
4.16 for this case having a Reynolds number of 1000. The recurring condition algorithm 
was used. As in the previous two time-dependent examples, the results lead directly to 
the steady solution after a sufficient number of time steps. This steady solution correlates 
closely with the results of Ghia et al (1982), as presented in Figure 4.15. 

It should be noted that Tosaka and Kakuda (1987) have run the transient driven cavity 
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at Re = 10,000. However, their results show signs of instability even at relatively small times, 
and are compared to the steady solution of Ghia et al which also is not correct at this 
much higher Reynolds number. A valid solution in this Re range would necessitate the use 
of an extremely refined mesh, far beyond that employed by Tosaka and Kakuda or Ghia 
et al. 

4.2.5. 7 Burgers Flow 

The classic uniaxial linear Burgers problem provides an excellent test of the convective 
thermoviscous formulations. The incompressible fluid flows in the x-direction with uniform 
velocity U. Meanwhile, the y-component of the velocity and temperature are specified as 
U 0 and T 0 , respectively, at inlet. Both are zero at the outlet. The length of the flow field 
is L. The analytical solution (Schlicting, 1955) is 

V„ = $U D 

T = < :T 0 

where 

? = {l -exp |fli, ( j- - l)] }/ {1- exp[-R L ]} 

with Rl = UL. 

The boundary element model employs eighteen quadratic surface elements encompass- 
ing the rectangular domain. The elements are graded, providing a very fine discretization 
near the exit, where V v and T vary substantially for large Rl. Results are shown in Figure 
4.17 for the thermal problem and in Figure 4.18 for the viscous problem. Excellent cor- 
relation with the analytical solution is obtained in both instances for this boundary-only 
analysis, even for the highly convective case of Rl = 1000. The portion of the flow field 
just ahead of the outlet is examined more closely in Figure 4.19. The convective Oseen 
solution obviously produces a precise solution. This problem can also be solved by utilizing 
the Stokes kernels and volume cells. As seen in Figure 4.19, this latter approach is not 
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quite as accurate. It should be noted that traditionally finite difference and finite element 
methods have a difficult time dealing with the convective terms present in this problem. 
Generally, ad hoc upwinding techniques must be introduced to produce stable, accurate 
solutions. On the other hand, with the convective boundary element approach the kernel 
functions contain an analytical form of upwinding. As a result, very precise BEM results 

can be obtained. 

4.2. 5.8 Flow Over a Cylinder 

As the final fluids example, the oft-studied case of incompressible flow over a circular 
cylinder is considered. In this problem, both the steady convective and non-convective 
formulations are utilized in the same analysis. The boundary element model is displayed 
in Figure 4.20. In the inner region, the Stokes kernels are employed along with a complete 
volume discretization. The outer region uses the Oseen kernels with a boundary-only 
formulation. The small non-linear contributions that would be present in the outer region 
away from the cylinder are ignored. The steady-state velocity vector plot at Re = 40 
' is shown in Figure 4.21. The recirculating zone, behind the cylinder, is clearly visible. 
Additionally, the resulting drag coefficient {C D ) of 1.8 obtained from the BE analysis is 
within the band of experimental scatter as presented by Panton (1984) for the circular 

cylinder. 

Interestingly, a completely linear Oseen analysis, which ignores all nonlinear convective 
terms in both regions, produces a very similar solution, except in the vicinity of the 
cylinder. Vector plots from the nonlinear analysis and the boundary-only linear Oseen 
analysis are superimposed in Figure 4.22. Although it is difficult to distinguish between the 
two analyses in that plot, both produce a recirculatory zone behind the cylinder. The linear 
solution, in general, overstates the velocities and velocity gradients in the neighborhood 
of the cylinder. Consequently, a drag coefficient of 3.4 is calculated, which is much higher 
than that found experimentally. This trend, of overpredicting the experimental drag, 
continues even to much higher Reynolds numbers as shown in Figure 4.23. Qualitatively, 
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however, the behavior of the BEM Oseen solution is consistent with the experimental 
curve for Reynolds Numbers up to 100,000. A much improved solution can be obtained by 
introducing a row of cells encompassing the cylinder. The full nonlinear problem is solved 
within this inner region, while the exterior remains a linear Oseen region. Figure 4.24 
illustrates a typical mesh, along with the resulting velocity vectors. As Reynolds number 
is increased, the significant nonlinear effects concentrate nearer to the cylinder, so that the 
thickness of the inner region may be reduced. Figure 4.23 also displays the drag obtained 
by utilizing just a single row of cells. Results are quite encouraging. Further improvement 
is possible by simply adding a second row of cells or by introducing a higher functional 
variation within each cell. The latter approach will be pursued in the coming year. 

4.3 Compressible Thermoviscous Flow 

4.3.1 Introduction 

Several of the previous examples have demonstrated the potential of the convective 
incompressible boundary integral formulation for flows in the high Reynolds number range. 
However, more generally, at very high speeds, compressibility of the fluid must also be 
considered. In particular, shock-related phenomena are not present in the incompressible 
formulations and kernel functions. To correct this deficiency, a compressible thermoviscous 
integral formulation is presented in this section. It should be noted that, while Oseen 
derived most of the fundamental solutions required for the incompressible case, no such 
similar solutions are available for compressibility. Consequently, considerable time and 
effort was required to derive these new approximate infinite space Green’s functions. A 
complete derivation of the compressible formulation was included in Dargush et al (1988) 
and will not be repeated here. This year only the main points will be highlighted, and the 
explicit form of the new kernel functions (Shi, 1990) will be provided. 

4.3.2 Governing Equations 

The conservation laws of mass, momentum and energy for a compressible thermovis- 
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cous fluid can be written in the following form 


dvi Dp , j n 

- f ipr Di + *=° 


(A + /x) 


d 2 v 


dxjdxi 

. a 2 e 


■} , 32 *>< 
— +M 


^- p ^i + /j=0 

dxjdxj dxi Dt 


DO dvi 

- pCtTZ - P— + 4> = 0 


(4.42a) 

(4.426) 

(4.42c) 


dxjdxj r ^ Dt dxi 

where <f> is a mass source and A is a second viscosity coefficient. All other quantities are 
defined in Section 4.2.2. Reference values for each of the primary variables are introduced 
in an effort to produce a linearized differential operator. Thus, let 


Vi = U, + u,- 

(4.43a) 

P = Vo + P 

(4.436) 

6 = 6 a + 6 

(4.43c) 

P = Po + P, 

(4.43d) 


in which and p Q are constant reference values, and Ui,p, 6 and p are the perturba- 

tions. Plugging these definitions into (4.42) produces, after some manipulation, 


dui D 0 p , 

-"•SS D, + *=° 

(4.44a) 

,, , \ d2u i , .. 32 “< d P . D 0 Ui , t , _ n 

dxjdxi +fi dxjdxj dxi P ° Dt + 

(4.446) 

, d 2 e d 0 s „ n 

k dx.dx , poCi Dt +rp ~ 

(4.44c) 


where and V' are now modified body mass sources, forces, and heat sources. Also, in 


(4.44), 


Do 

Dt 



(4.45) 
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4.3.3 Integral Representations 

4.3.3. 1 Steady Convective 

The formal appearance of the governing integral equations for steady compressible 
thermoviscous flow is very similar to that provided in Section 4.2.3.3. Specifically, let 

= J K,t. - fW dS + l l c M 

where once again 

«a = {«1 «2 0 } 

ta = {tl *2 q} 

f a = {/; fL va 

The major difference is, of course, in the kernel functions G ^ and • Actually, only the 
GVj and F? portion of these kernels changes with compressibility. The vortical component 
remains the same, however, there is now a dilatational component to the flow field that is 
absent in the incompressible case. 

This dilatational response is shown in Figures 4.25 for point forces in a compressible 
viscous fluid moving with a uniform reference velocity U x . The component <?u is plotted 
in the subsonic range in Figure 4.25a. The response increases as the magnitude of U\ is 
elevated. Notice also that G n is singular, with a sign reversal, at the load point. When U x 
equals the sonic velocity, there is a sudden dramatic change in the behavior of <?u, as can 
be seen in Figure 4.25b. Above the speed of sound, the response decreases and becomes 
more localized with increasing velocity. This is displayed in Figure 4.25c. On the other 
hand, G 22 which is a response perpendicular to the free stream, shows no discontinuity 
throughout the entire range of Ui. However, the response does peak at the sonic velocity. 
Evidently, these kernels do capture the nature of shock and, consequently, will be quite 
useful for the analysis of compressible thermoviscous flow. The complete kernel G^ derived 
by Shi (1990) is detailed in Appendix B.5. Implementation is planned for 1990. 
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FIGURES 4.2a,b 


CONVECTIVE THERMOV I 5C0U5 KERNELS 
C=l, u =1, U2=0. 0 , Y2=0. 0 










FIGURE 4.4a - CONVERGING CHANNEL 
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FIGURE 4.4b - CONVERGING CHANNEL 
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FIGURE 4.14 - DRIVEN CAVITY FI£W 
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DRIVEN CAVITY - FOUR REGION MODEL 
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HORIZONTAL VELOCITY at X=0.5 






FIGURE 4.16 - TRANSIENT DRIVEN CAVITY FLOW 
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FIGURE 4.21 
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FIGURE 4.25a,b 
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FIGURE 4.25c 
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5. FLUID-STRUCTURE INTERACTION 

5.1 Introduction 

In the previous two sections, boundary element formulations have been developed sep- 
arately for a thermoelastic structural component and for a thermoviscous fluid. However, 
the ultimate goal of this ongoing grant is to develop a single computer program to deter- 
mine the temperatures, deformation and stresses of a component exposed to a hot gas flow 
path, without the need for experimentally determined ambient fluid temperatures and film 
coefficients. While further work is still required for the fluid phase, sufficient progress has 
been made to demonstrate the utility of the overall concept. Consequently, in this section, 
problems of fluid-structure interaction will be examined. 

5.2 Formulation 

The Geometric Modeling Region (GMR) provides the vehicle for achieving interaction 
between the solid and fluid. Recall that in Section 4 different fluid formulations were 
employed in different GMRs. Now, some of the regions will use the thermoelastic solid 
boundary element model, while others utilize one of the thermoviscous fluid formulations. 
Compatibility must be enforced across all GMR interfaces, no matter which model is used 
for adjoining regions. 

For demonstration purposes, consider the problem of flow past a blade as sketched in 
Figure 5.1. The blade itself is labeled GMR1, and is modeled as a thermoelastic solid. 
A boundary mesh is all that is required for this structure. Surrounding the blade is a 
thin layer of cells. This is a nonlinear thermoviscous fluid region, named GMR2, which 
uses the Stokes formulation of Section 4.2.3.I. In this region, the complete Navier-Stokes 
equations are solved. Finally, the outer region GMR3 employs the convective Oseen kernels 
discussed in Section 4.2.3.3. Since no cells are present, the nonlinear volume and surface 
integrals in equation (4.13) are ignored. Thus, an approximation is introduced. However, as 
mentioned previously, outside of the boundary layer and wake these nonlinear contributions 
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are negligible. 

The interface between GMR2 and GMR3 poses no particular problem. Total velocity 
and temperature from both regions are equated at each interface node, while the tractions 
and flux must be equal in magnitude but of opposite direction. The latter conditions for 
the compatibility of traction and flux are also true for the solid-fluid interface between 
GMR1 and GMR2. Total temperature must, of course, be equal on this interface as well. 
However, the solid integral formulations of Section 3 are written in terms of displacement, 
while those for fluids use velocity. Consequently, a change in variable must be introduced 
to ensure complete interface compatibility. For that purpose, consider the following matrix 
form of the integral equation for a thermoviscous fluid. 

[•; il'M-P* il'M-l* £]'(?)*(*}• " 

The contributions from nonlinearities and past time steps are all contained in Re, as are 
any terms associated with the translation from perturbed velocity to total velocity 
Meanwhile, a similar expression written for a thermoelastic solid becomes 

[•; ;]'(•}- 1£ il’ftHR in?MS>’ " 

where u, is the total displacement. This must be rewritten in terms of total velocity v<, 

where 

dtn 


Vi = 


dr 


( 5 . 3 ) 


After invoking properties of the convolution integrals that are present in the original inte- 
gral equation (3.2), the appropriate representation for the solid can be written 

[? £]'{?}-[* M 

in which d i 3 ,d ej and P 0j are now modified kernel functions and k 0 is the corresponding 

right-hand-side contribution. However, at this point, the fluid formulation (5.1) and the 
solid formulation (5.4) are completely compatible, and are in an ideal form to solve quite 
general interaction problems. 
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5.3 Numerical Implementation 

The boundary element code, GPBEST, was generalized so that any combination of 
solid and fluid regions could be accommodated. Also, the modified thermoelastic kernels 
of equation (5.4) were implemented. The entire GPBEST input is free format and keyword 
driven. Output is provided on a region-by-region basis, and thus contains only informa- 
tion pertinent to the region type. Displacements, temperatures, stresses and strains are 
detailed for solid GMRs, while velocities, temperatures, stresses, pressures, strain rates 
and vorticities are output for fluid regions. In all cases, a complete PATRAN interface is 
available, so that any quantities can be plotted. 

5.4 Numerical Examples 

5.4.1 Introduction 

In this subsection a couple of examples will be presented to highlight the attractiveness 
of the present coupled boundary element approach. Flow past a thick-walled cylinder and 
an airfoil are considered. Both steady and transient conditions are examined, and a number 
of additional features of the GP-BEST implementation are explored. 

5.4.2 Steady Response of a T hick Cylinder 

For the first example, a thick-walled stainless steel cylinder rests under plane strain 
conditions in a stream of hot gas. The cylinder has an outer diameter of 1.0 in. and a 
thickness of 0.125 in. The inner surface of the cylinder is maintained at a temperature of 
0°F, while the gas temperature in the free stream is 1000°F. The following thermoelastic 
properties are assumed for the solid cylinder 

E = 29. X 10 6 pat, v = 0.30 

a = 9.6 x 10 -6 in./in °F 
k = 6.48 in.lb./sec.in.°F 

p = 7.34 x 10 -4 lb.sec. 2 /in. 4 c t = 3.83 x 10 5 in.lb.in./lb.sec. 2o F. 
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Additionally, the thermoviscous properties of the hot gas are taken as 
(i = 5.30 x 10~ 9 lb.sec/in. 2 

k = 7.28 x 10 -3 in.lb./sec.in. c> F 

p = 3.69 x 10 -8 lb.sec. 2 /in. 4 c p = 9.49 x 10 6 in.lb.m./lb.8ec. 2o F. 

Fluid velocities of 144 in. /sec., 1440 in. /sec. and 14400 in. /sec., corresponding to Reynolds 

Numbers of 10 3 , 10 4 and 10 5 , are examined. In all cases, the hot gas flows from left to right, 

and only the steady response is considered. 

At Re = 1000, the maximum temperature in the cylinder is only 98°F, and the peak 
compressive axial stress is 36 ksi. However, when the fluid velocity is increased to attain 
an Re = 10,000 a much more significant response is obtained. The temperature contours 
are shown in Figure 5.2a, the deformed shape is depicted in Figure 5.2b, and Figure 
5.2c illustrates the axial stress distribution. It should be noted that in Figure 5.2b the 
deformation has been scaled by a factor of 100. The effects of convection are quite evident 
in all three diagrams. With Reynolds number increased to 100,000 these effects .become 
even more pronounced, as seen in Figures 5.3. Now the peak metal temperature has 
reached 918°F. 

5.4.3 Airfoil Exposed to Hot Gas Flowpat h 

In this final example, an NACA0018 airfoil with an internal cooling passage is exposed 
to the flow of a hot gas. The boundary element model for the airfoil is shown in Fig- 
ure 5.4. Each dash represents an individual quadratic surface element. Throughout this 
problem, the outer gaseous region is modeled as a linear steady convective domain. Thus, 
a boundary-only exterior GMR is employed for the fluid. The hot gas at 1000°F flows 
from left to right, while the inner surface of the airfoil is maintained at 200°F . Material 
properties from the previous example are once again used to characterize both the solid 
and fluid. 

For the first set of investigations, the behavior of the airfoil is determined under steady- 
state conditions. Figure 5.5a displays the deformed shape at a Reynolds number of 1000 
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(based upon chord length). The solid line represents the final deformed shape, except 
that displacements have been scaled by a factor of twenty-five. Meanwhile, Figures 5.5b 
and c present the profiles of temperature and axial stress, respectively, along the upper 
surface of the airfoil. At this relatively slow speed flow, the airfoil is only effected near 
its leading edge. More significant response is shown in Figures 5.6a-c for Re = 10,000 and 
Figures 5.7a-c for Re = 100,000. In the latter case, the temperature at the stagnation point 
is nearly that of the free stream. All three cases considered so far have assumed an angle 
of attack of 0° with respect to the x-axis. Consequently, the response of the upper and 
lower surfaces is identical. Next, the angle of attack (a) is modified to 5° and 10°. Results 
for these cases are shown in Figures 5.8 and 5.9, respectively. Considerable asymmetry 
between upper and lower surfaces is now evident, although peak values of temperature and 
stress are essentially unaffected. 

Thermal barrier coatings are often employed to reduce the metal temperatures and 
stresses in hot section components. The benefit of such coatings can easily be evaluated 
with the present boundary element formulation. Consider, for example, a coating material 
with thermal conductivity Jfc = 0.50 in.lb./sec.in.°F sprayed to a thickness of .0095in. This is 
equivalent to an interfacial thermal resistance of .021 sec.in°F/in.lb., which can be specified 
on the fluid-to-solid GMR interface. Results are displayed in Figure 5.10 for Re = 100,000 
at a = 10°. Peak airfoil temperature is reduced from 976°F to 738°F by introducing this 
particular thermal barrier coating. 

Finally, it is of considerable interest to examine the transient response of the airfoil. 
At time zero, the airfoil is in thermal equilibrium at a temperature of 200°F. Suddenly, 
it is subjected to the hot gas stream with Re = 100,000 and a = 10°. The response of the 
upper surface at 1 msec., 2msec., 5 msec., and 10 msec, is shown in Figures 5.11-5.14. 
For this transient case, the peak stress occurs slightly offset from the tip of the airfoil. 
Additionally, the stress a vy reaches a maximum at approximately 2 msec., while a IX and 
the temperature continue to climb to their steady-state values. This is true of the axial 
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stress only because of the assumption of plane strain. In a full three-dimensional analysis, 
a „ would also have a higher peak during transient state. 


81 



FIGURE 5.1 - FLUID-STRUCTURE INTERACTION CONCEPTUAL MODEL 
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AIRFOIL - BOUNDARY ELEMENT MODEL 
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figure 5.9a - AIRFOIL (STEADY; Re = 100,000; 
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FIGURE 5.9b-e 
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AIRFOIL (STEADY: Re = 100,000; a = 10°) 
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figure 5.10 - AIRFOIL WITH COATING (STEADY; Re = 100,000; a = 10 










6. BEM FOR RELATED PHYSICAL PHENOMENA 

During the course of the investigation of the hot fluid-structure problem, a number of 
related technologies have been opened to analysis by the boundary element method. In this 
section, several of these potential applications are discussed. Most of the advancements 
depend upon the development of new fundamental solutions. For each case, a systematic 
procedure can be applied to obtain the required fundamental solution. This same procedure 
was developed and refined during the derivation of all of the kernel functions presented in 

Section 3 and 4. 

Perhaps the most interesting of these applications involve either moving sources or 
moving media. An example of the former kind is the determination of residual stresses in 
welds. As part of the NASA/HOST program, the boundary element code BEST3D was 
developed for the inelastic analysis of structures. Included in that code are a number of 
elastoplastic and viscoplastic material models that would be suitable for the weld problem. 
However, the temperature in the weld and adjoining structure is not known a priori, and 
a transient heat conduction analysis is required which accounts for the speed of the weld. 
The desired integral formulation for this thermal analysis is quite similar to that discussed 
for convective flow in Section 4. In addition, the fundamental solution that is needed for 
moving heat sources has already been derived as part of the present work. The other 
major advancement in boundary element technology that is required to solve the weld 
problem involves the development of more sophisticated nonlinear solution algorithms. It is 
envisioned that the modified Newton-Raphson schemes, employed for thermoviscous fluids, 
will provide the basis for that development. It should be noted that similar problems, such 
as frictional heating, grinding, and machining could also be studied utilizing the moving 
heat source approach. 

The hot viscous fluid formulations presented in Section 4 are quite general, and conse- 
quently, applicable to a wide range of physical processes. For example, the incompressible 
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integral equations could be used to solve the flow problem in injection molds, or the con- 
vective formulations could be applied to investigate the cooling of electronic components. 
Furthermore, some relatively minor extensions would provide significant beneflts. The in- 
clusion of a buoyancy term based upon the Boussinesq approximation, would permit the 
examination of the thermally-induced flow in lakes or the slow heating of a room. The 
addition of an extra equation involving the concentration of a diffusing substance provides 
the opportunity to investigate the spread of pollutants in a convective environment. 

As mentioned previously, once the techniques for obtaining fundamental solutions have 
been mastered, a wide range of physical phenomena can be analysed via boundary element 
approach. Recent work by Kaynia and Banerjee (1990) has focused on the development 
of fundamental solutions for dynamic poroelasticity. These solutions will be utilised in 
a BEM (Chen, 1991) for the analysis of soil-structure interaction under seismic loading. 
The analogous problem of dynamic thermoelasticity, which includes the important case of 
thermal shock, can also be solved with the same formulation. 

The coupling approach discussed in Section 6 can be used not only to solve the ther- 
moviscous fluid-structure problem, but also to investigate flutter. In this case, frequency- 
dependent formulation solutions are required. The infinite space solution for periodic 
elastodynamics of solids is well-known (Banerjee and Butterfield, 1981), while that for a 
linearized Oseen fluid could be derived. The frequency domain BEM analysis would be an 
extension of the work done for the NASA/HOST program and contained in BEST3D. 

There currently exists no satisfactory numerical nor analytical techniques to effectively 
deal with all of the physical phenomena mentioned in the preceding paragraphs. However, 
as an indirect result of the present hot fluid-structure grant, boundary element formulations 
and implementations are now possible for each case. 
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7. SUMMARY OF PROGRESS 


During this past year, major progress was made on two fronts. The first of these 
relates to the development of convective formulations for thermoviscous fluids. New kernels 
were derived in explicit form for both the incompressible and compressible cases. These 
kernels contain more of the physics of high speed flow than do the stationary kernels that 
were used previously. Consequently, high Reynolds number solutions are now possible, 
as demonstrated in several of the incompressible flow examples included in Section 4. 
Approximate solutions are obtained via either a linearized boundary-only model or by 
including volume cells just in the significantly nonlinear portions of the flow field. In 
order to obtain valid solutions at high Re , the numerical surface and volume integration 
algorithms within GP-BEST were completely revamped to handle the convective kernels. 
Meanwhile, the compressible formulation still requires some further work. 

The other major accomplishment concerns the implementation of a hot fluid-structure 
interaction capability. This capability was also added to GP-BEST in a very general man- 
ner so that considerable flexibility exists in terms of problem geometry, material properties, 
and boundary condition specification. Of particular significance is the fact that a single 
analysis code can now be used to analyze structures problems, fluids problems, or the 
complete fluid-structure interaction problem. As a result, transient thermal stresses in 
a hot section component can be obtained, at least approximately, without the need for 
experimentally determined convection coefficients and ambient temperatures. A couple of 
examples of the fluid-structure capability were provided in Section 5. Further enhancement 
of the BEM fluid formulation is still needed to improve these approximate solutions. 

A number of other important tasks were completed during this past year. Numerical 
integration was improved for the thermoelastic solid formulation. As a result, steady- 
state problems typically produce five-digit accuracy. The transient incompressible fluid 
formulation of Section 4. 2. 3. 2 was finalized. An algorithm was developed for boundary 
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pressure, vorticity, dilatation and stress in fluids. Routines for lift and drag calculation were 
added. A comprehensive PATRAN interface was written to permit the graphical display 
of all results. Lastly, a complete set of more than three dozen verification problems was 
created to test various facets of the GP-BEST code related to fluid-structure interaction. 
This ensures the maintenance of a reliable code, even during the ongoing development 
process. 


103 



8. WORKPLAN FOR THE NEXT YEAR 

Despite the significant progress that was made for high Reynolds Number flows in 
1989, the utility of the present program is still limited primarily by the ability to properly 
model the fluid which surrounds the hot section structural component of interest. The 
interaction facility, outlined in Section 5, is very general. Consequently, any number of 
boundary element solid models could be easily slotted into the current code. For example, 
the structure could be manufactured from a directionally-solidified material, or viscoplastic 
effects could be included. (These formulations have already been developed during the 
NASA/HOST program, although an extension to include transient thermal loading would 
be required.) Thus, most of the work that remains relates to fluid flow at high velocity. 

The workplan for the period from March 1990 to March 1991 is divided into two major 
tasks as outlined below. Task I can be completed with funding comparable to that provided 
during the present year. Task II requires an additional level of support. 

Task I 

• Implementation of the compressible convective formulation of Section 4.3. 

• Introduction of higher order cells for convective fluid regions, including full quadratic 
and/or quartic variation of velocity and temperature. 

• Development of a thin-GMR approximation for the boundary layer in high speed flow. 

• Semi-analytical integration of the convective kernels for singular and non-singular in- 
tegration in order to reduce the need for heavy numerical subsegmentation. 

• Development of a more efficient iterative solution algorithm for high speed flow. 
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Task II 


• Derivation of the required solid and fluid kernels for three-dimensional problems. 

• Implementation of the three-dimensional formulation within GP-BEST. The formula- 
tion will consist of a transient thermoelastic solid and an incompressible thermoviscous 
fluid using a steady convective fundamental solution. 

• Initial verification testing of the three-dimensional capability. 


t 
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APPENDIX B.l - 


Kernels for Thermoelasticity 


This appendix contains the detailed presentations of all the kernel functions utilized in 
the formulations contained in Section 3. Two-dimensional (plane strain) kernels are pro- 
vided, based upon continuous source and force fundamental solutions. For time-dependent 
uncoupled quasistatic thermoelasticity the following relationships must be used to deter- 
mine the proper form of the functions required in the boundary element discretization. 
That is, 

G n a0 {X - 0 = G a0 (X - e, nAt) for n - 1 

- £) = G a p{X - £, nAt) — G a pX - £, (n - l)At) for n > 1, 
with similar expressions holding for all the remaining kernels. In the specification of these 
kernels below, the arguments (X - £, t) are assumed. The indices 
vary from 1 to d 

a, p vary from 1 to (d + 1) 

6 equals d + 1 

where d is the dimensionality of the problem. Additionally, 

Xi coordinates of integration point 

& coordinates of field point 

y i = Zi-£i r' 2 =y i y i . 

For the displacement kernel, 


Gie = 0 

[(*)*«] 

G ** - h 0 ) 
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whereas, for the traction kernel, 


Fij ~ 4*r (1 


1 _ ^ SjjVkrik + Viny 'j ^ ^ 


+ 


(Sffl) (,-*,)] 


F„= 0 

f « = s(xfs)l(^) /sW - (n ' )/,w l 
*• - h K 2 ^) AW] • 


In the above, 



For the interior stress kernels, 


where 


BGjj 

96 


Epij = 

Dm = 


l l 
8ar n(l — v) 


[( 


2fiU fa 

aG *+ M ( 

<9G* 

+ 8 M 

— P&ijGpg 

1 - 2v 1 


< 9£y 

96 / 

s, 

dFpi 

.*7r + M 

(BFgj 

+ 8 j» 

- P&ijFpe 

1 — 2t/ 

J 96 

\ 9 6 

96 / 

2y*ykVk 

r 3 

1 

"t ** 
£ 

1 

1 On 

f") 

+ /Mi' 

| (3 - 4*)] 
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dFij _ 1 1 [_ / 4y,yyt/fcy t ni _ y,y 3 n k _ 6 3k yiyin t 

dT k ~ 4 tit 2 (l — i/) [_ \ r4 r2 r2 

_ Sji^n, j Aw _ (few _ ^ + fe> - /,(„) 

+ - fe<) Ad)] 

2^i = J_ f _£_\ f /' few ') { 2 S, _ ,-.*/n _ (»5£ + E2i + few') (t,)' . 
d£ fc 4*r V* + 2 m/ [V r 3 / X V r r r ) J 

AM = 2 

AM = 1-21/ 

AM = 1-21/ 
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APPENDIX B.2 - 

Kernels for Steady Incompressible Thermoviscous Flow 




Fa = 


1 [ 2y,yjVfcrtfc l 

2*r[ r 3 J 


a<?,y _ 1 f gyfcy» | ftfcVj _ ftjVfc 

<9xfc 47r^ir r r r 




Fee — 


l \ ykKk 
2irr [ r J 


dG$e _ 1 Vk 

dxk 2nkr r m 


Vi = ~ ti 


r 2 = yiVi 


2 ViViVk 
r 3 
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APPENDIX B.3 - 

Kernels for Unsteady Incompressible Viscous Flow 


1 Tv<V,, , ,, r 


F.M-XA = £{*?*(•><■.> - '-''“I + *?<’■<’> - m) + 




where 


_ 2y»yg'fc» fc^ 23i ( r? ) _ e -n J /<} 


dG <i ( Z-X,t) = — 


dx k 


Airfir 


fa&O.M} + ^M*» - 

r r " 


-h&pSLiuM 


y> = &•■*» r2 — ww 

»?=(5fl7T C = ^f> 

«i(»?) = - e ~ V ^ ) 

Ei{z) = l°° ^dxx. 


Then, 


G£(e - X) = GiM - X, nAr) fo r 

G?j(t ~X)= G,;(£ - X, nAr) - G i} {{ -X,(n- l)Ar) for 


with similar relationships for — X) and 0xi j (£ •X’) 


*-sx{v)} 


n = 1 
n > 1 
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APPENDIX B.4 - 

Kernels for Steady Convective Incompressible Viscous Flow 


Gij - 



d<t> c (ua a* c 

dx 3 U\u)dxi U 



Gpj - 2* 




where 


yi = X\ - i ^ ~ 

c = n IP = UiUi 

p 

P = U k y k /2c 
a = Ur /2c 


d<j> 

dxi 


4> = — in(a) — e ^ K 0 {a ) 
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APPENDIX B.5 - 


Kernels for Steady Convective Compressible Viscous Flow 


where 


G i} = ^ 6 t] e- 0 K o (a) + ~(U,y 3 + U jVi - 6 i} U k y k ) e ~^ K^a) 


H(c, -U) c 1 
2np 0 V U 2 R 2 [ 


UiVi + Uj-ffi - 6 i} U k y k + U -^U iU} 


?2 _ ,2 


H(c. 


yi = x« - ii 

£ 

5 

II 

C = /i/p 

U 2 = UiUi 

0 = u k y k /2c 

a = Ur/2c 

C 2 . = Po/Po 

V 2 = c 2 - U 2 

■♦W 


f 1 for U < c, 
~ \ 0 for U > c. 

subsonic 

sonic and supersonic 
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